/******************************************************************************
 * 
 * Announce: CSharpKit, Basic algorithms, components and definitions.
 *           Copyright (C) ShenYongchen.
 *           All rights reserved.
 *   Author: 申永辰.郑州 (shenyczz@163.com)
 *  WebSite: http://gitee.com/shenyczz/csharp_csk
 *
 * THIS CODE IS LICENSED UNDER THE MIT LICENSE (MIT).
 * THIS CODE IS PROVIDED *AS IS* WITHOUT WARRANTY OF 
 * ANY KIND, EITHER EXPRESS OR IMPLIED, INCLUDING ANY
 * IMPLIED WARRANTIES OF FITNESS FOR A PARTICULAR
 * PURPOSE, MERCHANTABILITY, OR NON-INFRINGEMENT.
 * 
******************************************************************************/

using System;
using System.Numerics;

namespace CSharpKit.Numerics
{
    // References:
    // [1] https://github.com/scipy/scipy/blob/master/scipy/special/amos_wrappers.c
    public static class Amos
    {
        #region AiryAi

        // Returns Ai(z)
        public static Complex Cairy(Complex z)
        {
            int id = 0;
            int kode = 1;
            int nz = 0;
            int ierr = 0;

            double air = double.NaN;
            double aii = double.NaN;

            AmosHelper.zairy(z.Real, z.Imaginary, id, kode, ref air, ref aii, ref nz, ref ierr);
            return new Complex(air, aii);
        }

        // Returns Exp(zta) * Ai(z) where zta = (2/3) * z * Sqrt(z)
        public static Complex ScaledCairy(Complex z)
        {
            int id = 0;
            int kode = 2;
            int nz = 0;
            int ierr = 0;

            double air = double.NaN;
            double aii = double.NaN;

            AmosHelper.zairy(z.Real, z.Imaginary, id, kode, ref air, ref aii, ref nz, ref ierr);
            return new Complex(air, aii);
        }

        // Returns Exp(zta) * Ai(z) where zta = (2/3) * z * Sqrt(z)
        public static double ScaledCairy(double z)
        {
            if (z < 0)
            {
                return double.NaN;
            }

            int id = 0;
            int kode = 2;
            int nz = 0;
            int ierr = 0;

            double air = double.NaN;
            double aii = double.NaN;

            AmosHelper.zairy(z, 0.0, id, kode, ref air, ref aii, ref nz, ref ierr);
            return air;
        }

        // Returns d/dz Ai(z)
        public static Complex CairyPrime(Complex z)
        {
            int id = 1;
            int kode = 1;
            int nz = 0;
            int ierr = 0;

            double air = double.NaN;
            double aii = double.NaN;

            AmosHelper.zairy(z.Real, z.Imaginary, id, kode, ref air, ref aii, ref nz, ref ierr);
            return new Complex(air, aii);
        }

        // Returns Exp(zta) * d/dz Ai(z) where zta = (2/3) * z * Sqrt(z)
        public static Complex ScaledCairyPrime(Complex z)
        {
            int id = 1;
            int kode = 2;
            int nz = 0;
            int ierr = 0;

            double air = double.NaN;
            double aii = double.NaN;

            AmosHelper.zairy(z.Real, z.Imaginary, id, kode, ref air, ref aii, ref nz, ref ierr);
            return new Complex(air, aii);
        }

        // Returns Exp(zta) * d/dz Ai(z) where zta = (2/3) * z * Sqrt(z)
        public static double ScaledCairyPrime(double z)
        {
            if (z < 0)
            {
                return double.NaN;
            }

            int id = 1;
            int kode = 2;
            int nz = 0;
            int ierr = 0;

            double air = double.NaN;
            double aii = double.NaN;

            AmosHelper.zairy(z, 0.0, id, kode, ref air, ref aii, ref nz, ref ierr);
            return air;
        }

        #endregion

        #region AiryBi

        // Returns Bi(z)
        public static Complex Cbiry(Complex z)
        {
            int id = 0;
            int kode = 1;
            int nz = 0;
            int ierr = 0;

            double bir = double.NaN;
            double bii = double.NaN;

            AmosHelper.zbiry(z.Real, z.Imaginary, id, kode, ref bir, ref bii, ref nz, ref ierr);
            return new Complex(bir, bii);
        }

        // Returns Exp(-axzta) * Bi(z) where zta = (2 / 3) * z * Sqrt(z) and axzta = Abs(zta.Real)
        public static Complex ScaledCbiry(Complex z)
        {
            int id = 0;
            int kode = 2;
            int nz = 0;
            int ierr = 0;

            double bir = double.NaN;
            double bii = double.NaN;

            AmosHelper.zbiry(z.Real, z.Imaginary, id, kode, ref bir, ref bii, ref nz, ref ierr);
            return new Complex(bir, bii);
        }

        // Returns d/dz Bi(z)
        public static Complex CbiryPrime(Complex z)
        {
            int id = 1;
            int kode = 1;
            int nz = 0;
            int ierr = 0;

            double bipr = double.NaN;
            double bipi = double.NaN;

            AmosHelper.zbiry(z.Real, z.Imaginary, id, kode, ref bipr, ref bipi, ref nz, ref ierr);
            return new Complex(bipr, bipi);
        }

        // Returns Exp(-axzta) * d/dz Bi(z) where zta = (2 / 3) * z * Sqrt(z) and axzta = Abs(zta.Real)
        public static Complex ScaledCbiryPrime(Complex z)
        {
            int id = 1;
            int kode = 2;
            int nz = 0;
            int ierr = 0;

            double bipr = double.NaN;
            double bipi = double.NaN;

            AmosHelper.zbiry(z.Real, z.Imaginary, id, kode, ref bipr, ref bipi, ref nz, ref ierr);
            return new Complex(bipr, bipi);
        }

        #endregion

        #region BesselJ

        // Return J(v, z)
        public static Complex Cbesj(double v, Complex z)
        {
            if (double.IsNaN(v) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary))
            {
                return new Complex(double.NaN, double.NaN);
            }

            int sign = 1;
            if (v < 0)
            {
                v = -v;
                sign = -1;
            }

            int n = 1;
            int kode = 1;
            int nz = 0;
            int ierr = 0;
            double[] cyjr = new double[n];
            double[] cyji = new double[n];
            for (int i = 0; i < n; i++)
            {
                cyjr[i] = double.NaN;
                cyji[i] = double.NaN;
            }

            AmosHelper.zbesj(z.Real, z.Imaginary, v, kode, n, cyjr, cyji, ref nz, ref ierr);
            Complex cyj = new Complex(cyjr[0], cyji[0]);

            if (ierr == 2)
            {
                //overflow
                cyj = ScaledCbesj(v, z);
                cyj = new Complex(cyj.Real * double.PositiveInfinity, cyj.Imaginary * double.PositiveInfinity);
            }

            if (sign == -1)
            {
                if (!ReflectJY(ref cyj, v))
                {
                    double[] cyyr = new double[n];
                    double[] cyyi = new double[n];
                    double[] cwrkr = new double[n];
                    double[] cwrki = new double[n];

                    for (int i = 0; i < n; i++)
                    {
                        cyyr[i] = double.NaN;
                        cyyi[i] = double.NaN;
                        cwrkr[i] = double.NaN;
                        cwrki[i] = double.NaN;
                    }

                    AmosHelper.zbesy(z.Real, z.Imaginary, v, kode, n, cyyr, cyyi, ref nz, cwrkr, cwrki, ref ierr);
                    Complex cyy = new Complex(cyyr[0], cyyi[0]);

                    cyj = RotateJY(cyj, cyy, v);
                }
            }
            return cyj;
        }

        // Return J(v, z)
        public static double Cbesj(double v, double z)
        {
            if (z < 0 && v != (int)v)
            {
                return double.NaN;
            }

            return Cbesj(v, new Complex(z, 0)).Real;
        }

        // Return Exp(-Abs(y)) * J(v, z) where y = z.Imaginary
        public static Complex ScaledCbesj(double v, Complex z)
        {
            if (double.IsNaN(v) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary))
            {
                return new Complex(double.NaN, double.NaN);
            }

            int sign = 1;
            if (v < 0)
            {
                v = -v;
                sign = -1;
            }

            int n = 1;
            int kode = 2;
            int nz = 0;
            int ierr = 0;

            double[] cyjr = new double[n];
            double[] cyji = new double[n];
            for (int i = 0; i < n; i++)
            {
                cyjr[i] = double.NaN;
                cyji[i] = double.NaN;
            }

            AmosHelper.zbesj(z.Real, z.Imaginary, v, kode, n, cyjr, cyji, ref nz, ref ierr);
            Complex cyj = new Complex(cyjr[0], cyji[0]);

            if (sign == -1)
            {
                if (!ReflectJY(ref cyj, v))
                {
                    double[] cyyr = new double[n];
                    double[] cyyi = new double[n];
                    double[] cworkr = new double[n];
                    double[] cworki = new double[n];
                    for (int i = 0; i < n; i++)
                    {
                        cyyr[i] = double.NaN;
                        cyyi[i] = double.NaN;
                        cworkr[i] = double.NaN;
                        cworki[i] = double.NaN;
                    }

                    AmosHelper.zbesy(z.Real, z.Imaginary, v, kode, n, cyyr, cyyi, ref nz, cworkr, cworki, ref ierr);
                    Complex cyy = new Complex(cyyr[0], cyyi[0]);

                    cyj = RotateJY(cyj, cyy, v);
                }
            }
            return cyj;
        }

        // Return Exp(-Abs(y)) * J(v, z) where y = z.Imaginary
        public static double ScaledCbesj(double v, double z)
        {
            if (z < 0 && v != (int)v)
            {
                return double.NaN;
            }

            return ScaledCbesj(v, new Complex(z, 0)).Real;
        }

        #endregion

        #region BesselY

        // Return Y(v, z)
        public static Complex Cbesy(double v, Complex z)
        {
            if (double.IsNaN(v) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary))
            {
                return new Complex(double.NaN, double.NaN);
            }

            int sign = 1;
            if (v < 0)
            {
                v = -v;
                sign = -1;
            }

            int n = 1;
            int kode = 1;
            int nz = 0;
            int ierr = 0;
            Complex cyy;

            if (z.Real == 0 && z.Imaginary == 0)
            {
                //overflow
                cyy = new Complex(double.NegativeInfinity, 0);
            }
            else
            {
                double[] cyyr = new double[n];
                double[] cyyi = new double[n];
                double[] cworkr = new double[n];
                double[] cworki = new double[n];
                for (int i = 0; i < n; i++)
                {
                    cyyr[i] = double.NaN;
                    cyyi[i] = double.NaN;
                    cworkr[i] = double.NaN;
                    cworki[i] = double.NaN;
                }

                AmosHelper.zbesy(z.Real, z.Imaginary, v, kode, n, cyyr, cyyi, ref nz, cworkr, cworki, ref ierr);
                cyy = new Complex(cyyr[0], cyyi[0]);

                if (ierr == 2)
                {
                    if (z.Real >= 0 && z.Imaginary == 0)
                    {
                        //overflow
                        cyy = new Complex(double.NegativeInfinity, 0);
                    }
                }
            }

            if (sign == -1)
            {
                if (!ReflectJY(ref cyy, v))
                {
                    double[] cyjr = new double[n];
                    double[] cyji = new double[n];
                    for (int i = 0; i < n; i++)
                    {
                        cyjr[i] = double.NaN;
                        cyji[i] = double.NaN;
                    }

                    AmosHelper.zbesj(z.Real, z.Imaginary, v, kode, n, cyjr, cyji, ref nz, ref ierr);
                    Complex cyj = new Complex(cyjr[0], cyji[0]);

                    cyy = RotateJY(cyy, cyj, -v);
                }
            }
            return cyy;
        }

        // Return Y(v, z)
        public static double Cbesy(double v, double x)
        {
            if (x < 0.0)
            {
                return double.NaN;
            }

            Complex z = new Complex(x, 0.0);
            return Cbesy(v, z).Real;
        }

        // Return Exp(-Abs(y)) * Y(v, z) where y = z.Imaginary
        public static Complex ScaledCbesy(double v, Complex z)
        {
            if (double.IsNaN(v) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary))
            {
                return new Complex(double.NaN, double.NaN);
            }

            int sign = 1;
            if (v < 0)
            {
                v = -v;
                sign = -1;
            }

            int n = 1;
            int kode = 2;
            int nz = 0;
            int ierr = 0;

            double[] cyyr = new double[n];
            double[] cyyi = new double[n];
            double[] cwrkr = new double[n];
            double[] cwrki = new double[n];
            for (int i = 0; i < n; i++)
            {
                cyyr[i] = double.NaN;
                cyyi[i] = double.NaN;
                cwrkr[i] = double.NaN;
                cwrki[i] = double.NaN;
            }

            AmosHelper.zbesy(z.Real, z.Imaginary, v, kode, n, cyyr, cyyi, ref nz, cwrkr, cwrki, ref ierr);
            Complex cyy = new Complex(cyyr[0], cyyi[0]);

            if (ierr == 2)
            {
                if (z.Real >= 0 && z.Imaginary == 0)
                {
                    //overflow
                    cyy = new Complex(double.PositiveInfinity, 0);
                }
            }

            if (sign == -1)
            {
                if (!ReflectJY(ref cyy, v))
                {
                    double[] cyjr = new double[n];
                    double[] cyji = new double[n];
                    for (int i = 0; i < n; i++)
                    {
                        cyjr[i] = double.NaN;
                        cyji[i] = double.NaN;
                    }

                    AmosHelper.zbesj(z.Real, z.Imaginary, v, kode, n, cyjr, cyji, ref nz, ref ierr);
                    Complex cyj = new Complex(cyjr[0], cyji[0]);
                    cyy = RotateJY(cyy, cyj, -v);
                }
            }
            return cyy;
        }

        // Return Exp(-Abs(y)) * Y(v, z) where y = z.Imaginary
        public static double ScaledCbesy(double v, double x)
        {
            if (x < 0)
            {
                return double.NaN;
            }

            return ScaledCbesy(v, new Complex(x, 0)).Real;
        }

        #endregion

        #region BesselI

        // Returns I(v, z)
        public static Complex Cbesi(double v, Complex z)
        {
            if (double.IsNaN(v) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary))
            {
                return new Complex(double.NaN, double.NaN);
            }

            int sign = 1;
            if (v < 0)
            {
                v = -v;
                sign = -1;
            }

            int n = 1;
            int kode = 1;
            int nz = 0;
            int ierr = 0;

            double[] cyir = new double[n];
            double[] cyii = new double[n];
            for (int i = 0; i < n; i++)
            {
                cyir[i] = double.NaN;
                cyii[i] = double.NaN;
            }

            AmosHelper.zbesi(z.Real, z.Imaginary, v, kode, n, cyir, cyii, ref nz, ref ierr);
            Complex cyi = new Complex(cyir[0], cyii[0]);

            if (ierr == 2)
            {
                //overflow
                if (z.Imaginary == 0 && (z.Real >= 0 || v == System.Math.Floor(v)))
                {
                    if (z.Real < 0 && v / 2 != System.Math.Floor(v / 2))
                        cyi = new Complex(double.NegativeInfinity, 0);
                    else
                        cyi = new Complex(double.PositiveInfinity, 0);
                }
                else
                {
                    cyi = ScaledCbesi(v * sign, z);
                    cyi = new Complex(cyi.Real * double.PositiveInfinity, cyi.Imaginary * double.PositiveInfinity);
                }
            }

            if (sign == -1)
            {
                if (!ReflectI(cyi, v))
                {
                    double[] cykr = new double[n];
                    double[] cyki = new double[n];
                    AmosHelper.zbesk(z.Real, z.Imaginary, v, kode, n, cykr, cyki, ref nz, ref ierr);
                    Complex cyk = new Complex(cykr[0], cyki[0]);

                    cyi = RotateI(cyi, cyk, v);
                }
            }

            return cyi;
        }

        // Return Exp(-Abs(x)) * I(v, z) where x = z.Real
        public static Complex ScaledCbesi(double v, Complex z)
        {
            if (double.IsNaN(v) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary))
            {
                return new Complex(double.NaN, double.NaN);
            }

            int sign = 1;
            if (v < 0)
            {
                v = -v;
                sign = -1;
            }

            int n = 1;
            int kode = 2;
            int nz = 0;
            int ierr = 0;

            double[] cyir = new double[n];
            double[] cyii = new double[n];
            for (int i = 0; i < n; i++)
            {
                cyir[i] = double.NaN;
                cyii[i] = double.NaN;
            }

            AmosHelper.zbesi(z.Real, z.Imaginary, v, kode, n, cyir, cyii, ref nz, ref ierr);
            Complex cyi = new Complex(cyir[0], cyii[0]);

            if (sign == -1)
            {
                if (!ReflectI(cyi, v))
                {
                    double[] cykr = new double[n];
                    double[] cyki = new double[n];
                    AmosHelper.zbesk(z.Real, z.Imaginary, v, kode, n, cykr, cyki, ref nz, ref ierr);
                    Complex cyk = new Complex(cykr[0], cyki[0]);

                    //adjust scaling to match zbesi
                    cyk = Rotate(cyk, -z.Imaginary / System.Math.PI);
                    if (z.Real > 0)
                    {
                        cyk = new Complex(cyk.Real * System.Math.Exp(-2 * z.Real), cyk.Imaginary * System.Math.Exp(-2 * z.Real));
                    }
                    //v -> -v
                    cyi = RotateI(cyi, cyk, v);
                }
            }

            return cyi;
        }

        // Return Exp(-Abs(x)) * I(v, z) where x = z.Real
        public static double ScaledCbesi(double v, double x)
        {
            if (v != System.Math.Floor(v) && x < 0)
            {
                return double.NaN;
            }

            return ScaledCbesi(v, new Complex(x, 0)).Real;
        }

        #endregion

        #region BesselK

        // Returns K(v, z)
        public static Complex Cbesk(double v, Complex z)
        {
            if (double.IsNaN(v) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary))
            {
                return new Complex(double.NaN, double.NaN);
            }
            if (v < 0)
            {
                //K_v == K_{-v} even for non-integer v
                v = -v;
            }

            int n = 1;
            int kode = 1;
            int nz = 0;
            int ierr = 0;

            double[] cykr = new double[n];
            double[] cyki = new double[n];
            for (int i = 0; i < n; i++)
            {
                cykr[i] = double.NaN;
                cyki[i] = double.NaN;
            }

            AmosHelper.zbesk(z.Real, z.Imaginary, v, kode, n, cykr, cyki, ref nz, ref ierr);
            Complex cyk = new Complex(cykr[0], cyki[0]);

            if (ierr == 1)
            {
                if (z.Real == 0.0 && z.Imaginary == 0.0)
                {
                    cyk = new Complex(double.PositiveInfinity, 0);
                }
            }
            else if (ierr == 2)
            {
                if (z.Real >= 0 && z.Imaginary == 0)
                {
                    //overflow
                    cyk = new Complex(double.PositiveInfinity, 0);
                }
            }

            return cyk;
        }

        // Returns K(v, z)
        public static double Cbesk(double v, double z)
        {
            if (z < 0)
            {
                return double.NaN;
            }
            else if (z == 0)
            {
                return double.PositiveInfinity;
            }
            else if (z > 710 * (1 + System.Math.Abs(v)))
            {
                // Underflow. See uniform expansion https://dlmf.nist.gov/10.41
                // This condition is not a strict bound (it can underflow earlier),
                // rather, we are here working around a restriction in AMOS.

                return 0;
            }
            else
            {
                Complex w = new Complex(z, 0);
                return Cbesk(v, w).Real;
            }
        }

        // returns Exp(z) * K(v, z)
        public static Complex ScaledCbesk(double v, Complex z)
        {
            if (double.IsNaN(v) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary))
            {
                return new Complex(double.NaN, double.NaN);
            }

            if (v < 0)
            {
                //K_v == K_{-v} even for non-integer v
                v = -v;
            }

            int n = 1;
            int kode = 2;
            int nz = 0;
            int ierr = 0;

            double[] cykr = new double[n];
            double[] cyki = new double[n];
            for (int i = 0; i < n; i++)
            {
                cykr[i] = double.NaN;
                cyki[i] = double.NaN;
            }

            AmosHelper.zbesk(z.Real, z.Imaginary, v, kode, n, cykr, cyki, ref nz, ref ierr);
            Complex cyk = new Complex(cykr[0], cyki[0]);

            if (ierr == 2)
            {
                if (z.Real >= 0 && z.Imaginary == 0)
                {
                    //overflow
                    cyk = new Complex(double.PositiveInfinity, 0);
                }
            }

            return cyk;
        }

        // returns Exp(z) * K(v, z)
        public static double ScaledCbesk(double v, double z)
        {
            if (z < 0)
            {
                return double.NaN;
            }
            else if (z == 0)
            {
                return double.PositiveInfinity;
            }
            else
            {
                Complex w = new Complex(z, 0);
                return ScaledCbesk(v, w).Real;
            }
        }

        #endregion

        #region HankelH1

        // Returns H1(v, z)
        public static Complex Cbesh1(double v, Complex z)
        {
            if (double.IsNaN(v) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary))
            {
                return new Complex(double.NaN, double.NaN); ;
            }

            int n = 1;
            int kode = 1;
            int m = 1;
            int nz = 0;
            int ierr = 0;

            double[] cyhr = new double[n];
            double[] cyhi = new double[n];
            for (int i = 0; i < n; i++)
            {
                cyhr[i] = double.NaN;
                cyhi[i] = double.NaN;
            }

            int sign = 1;
            if (v < 0)
            {
                v = -v;
                sign = -1;
            }

            AmosHelper.zbesh(z.Real, z.Imaginary, v, kode, m, n, cyhr, cyhi, ref nz, ref ierr);
            Complex cyh = new Complex(cyhr[0], cyhi[0]);

            if (sign == -1)
            {
                cyh = Rotate(cyh, v);
            }

            return cyh;
        }

        // Returns Exp(-z * j) * H1(n, z) where j = Sqrt(-1)
        public static Complex ScaledCbesh1(double v, Complex z)
        {
            if (double.IsNaN(v) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary))
            {
                return new Complex(double.NaN, double.NaN); ;
            }

            int n = 1;
            int kode = 2;
            int m = 1;
            int nz = 0;
            int ierr = 0;

            double[] cyhr = new double[n];
            double[] cyhi = new double[n];
            for (int i = 0; i < n; i++)
            {
                cyhr[i] = double.NaN;
                cyhi[i] = double.NaN;
            }

            int sign = 1;
            if (v < 0)
            {
                v = -v;
                sign = -1;
            }

            AmosHelper.zbesh(z.Real, z.Imaginary, v, kode, m, n, cyhr, cyhi, ref nz, ref ierr);
            Complex cyh = new Complex(cyhr[0], cyhi[0]);

            if (sign == -1)
            {
                cyh = Rotate(cyh, v);
            }

            return cyh;
        }

        #endregion

        #region HankelH2

        // Returns H2(v, z)
        public static Complex Cbesh2(double v, Complex z)
        {
            if (double.IsNaN(v) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary))
            {
                return new Complex(double.NaN, double.NaN);
            }

            if (v == 0 && z.Real == 0 && z.Imaginary == 0)
            {
                return new Complex(double.NaN, double.NaN); // ComplexInfinity
            }

            int n = 1;
            int kode = 1;
            int m = 2;
            int nz = 0;
            int ierr = 0;

            double[] cyhr = new double[n];
            double[] cyhi = new double[n];
            for (int i = 0; i < n; i++)
            {
                cyhr[i] = double.NaN;
                cyhi[i] = double.NaN;
            }

            int sign = 1;
            if (v < 0)
            {
                v = -v;
                sign = -1;
            }

            AmosHelper.zbesh(z.Real, z.Imaginary, v, kode, m, n, cyhr, cyhi, ref nz, ref ierr);
            Complex cyh = new Complex(cyhr[0], cyhi[0]);

            if (sign == -1)
            {
                cyh = Rotate(cyh, -v);
            }
            return cyh;
        }

        // Returns Exp(z * j) * H2(n, z) where j = Sqrt(-1)
        public static Complex ScaledCbesh2(double v, Complex z)
        {
            if (double.IsNaN(v) || double.IsNaN(z.Real) || double.IsNaN(z.Imaginary))
            {
                return new Complex(double.NaN, double.NaN);
            }

            if (v == 0 && z.Real == 0 && z.Imaginary == 0)
            {
                return new Complex(double.NaN, double.NaN); // ComplexInfinity
            }

            int n = 1;
            int kode = 2;
            int m = 2;
            int nz = 0;
            int ierr = 0;

            double[] cyhr = new double[n];
            double[] cyhi = new double[n];
            for (int i = 0; i < n; i++)
            {
                cyhr[i] = double.NaN;
                cyhi[i] = double.NaN;
            }

            int sign = 1;
            if (v < 0)
            {
                v = -v;
                sign = -1;
            }

            AmosHelper.zbesh(z.Real, z.Imaginary, v, kode, m, n, cyhr, cyhi, ref nz, ref ierr);
            Complex cyh = new Complex(cyhr[0], cyhi[0]);

            if (sign == -1)
            {
                cyh = Rotate(cyh, -v);
            }
            return cyh;
        }

        #endregion


        #region utilities

        /// <summary>
        /// return sin(pi * x)
        /// </summary>
        /// <param name="x"></param>
        /// <returns></returns>
        private static double SinPi(double x)
        {
            if (System.Math.Floor(x) == x && System.Math.Abs(x) < 1.0e14)
            {
                //Return 0 when at exact zero, as long as the floating point number is
                //small enough to distinguish integer points from other points.

                return 0;
            }
            return System.Math.Sin(System.Math.PI * x);
        }
        /// <summary>
        /// return cos(pi * x)
        /// </summary>
        /// <param name="x"></param>
        /// <returns></returns>
        private static double CosPi(double x)
        {
            if (System.Math.Floor(x + 0.5) == x + 0.5 && System.Math.Abs(x) < 1.0E14)
            {
                //Return 0 when at exact zero, as long as the floating point number is
                //small enough to distinguish integer points from other points.

                return 0;
            }
            return System.Math.Cos(System.Math.PI * x);
        }

        /// <summary>
        /// 旋转
        /// </summary>
        /// <param name="z"></param>
        /// <param name="v"></param>
        /// <returns></returns>
        private static Complex Rotate(Complex z, double v)
        {
            double c = CosPi(v);
            double s = SinPi(v);
            return new Complex(z.Real * c - z.Imaginary * s, z.Real * s + z.Imaginary * c);
        }

        private static Complex RotateJY(Complex j, Complex y, double v)
        {
            double c = CosPi(v);
            double s = SinPi(v);
            return new Complex(j.Real * c - y.Real * s, j.Imaginary * c - y.Imaginary * s);
        }

        private static bool ReflectJY(ref Complex jy, double v)
        {
            //NB: Y_v may be huge near negative integers -- so handle exact
            //     integers carefully

            if (v != System.Math.Floor(v))
            {
                return false;
            }

            int i = (int)(v - 16384.0 * System.Math.Floor(v / 16384.0));
            if (i % 2 == 1)
            {
                jy = new Complex(-jy.Real, -jy.Imaginary);
            }

            return true;
        }

        private static bool ReflectI(Complex ik, double v)
        {
            if (v != System.Math.Floor(v))
            {
                return false;
            }

            return true; //I is symmetric for integer v
        }

        private static Complex RotateI(Complex i, Complex k, double v)
        {
            double s = System.Math.Sin(v * System.Math.PI) * (2.0 / System.Math.PI);
            return new Complex(i.Real + s * k.Real, i.Imaginary + s * k.Imaginary);
        }

        #endregion

        //}}@@@
    }





}
